# -*- coding: utf-8 -*-
"""
Created on Thu Jul 28 08:18:36 2022

@author: Derek Rodriguez Pacheco, PhD 
IUSS Pavia
"""

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import subprocess
from scipy.optimize import minimize


## Model's parameters

betaf= 0.51
pinchX= 0.5
pinchY= 0.01

params = [betaf,pinchX,pinchY]


def HystModelFit(par):
    
    
    m= 950                           ## Mass in N*s2/m or kg
    g= 9.81                          ## Acceleration of gravity in m/s2
    
    adjfactors= np.array(pd.read_csv('Adjustment_Factors_Displacement_CM.txt', header= None, delim_whitespace= True)[0])
    indexesb= np.array(pd.read_csv('Indexes_Beginning.txt', header= None, delim_whitespace= True)[0])
    indexesf= np.array(pd.read_csv('Indexes_End.txt', header= None, delim_whitespace= True)[0])
    
    k1= 0.5
    k2= 0.01
    sigAct= 3.8
    epsSlip= 0
    epsBear= 0
    rBear= 1
    
    s1p= 2
    e1p= 6
    s2p= 4
    e2p= 25
    s3p= 15
    e3p= 20
    
    damage1= 0
    damage2= 0
    beta= 0
    
    s1n= s1p*31
    e1n= np.copy(e1p)
    s2n= s2p*31
    e2n= np.copy(e2p)
    s3n= s3p*31
    e3n= np.copy(e3p)
    
    
    width= 600
    height= 1053
    
    testnumber= [13,14,15,16,17,18,19,20,24,25,26,27,28,40,41]
    
    ##
    
    datacm= {}
    # acceltcm01= {}
    # acceltcm02= {}
    accelcm= {}
    timecm= {}
    timecmc= {}
    
    data= {}
    # accelda1= {}
    # accelro= {}
    accelbase= {}
    dispdacm= {}
    timearray= {}
    energyarraycm= {}
    
    timetotal= []
    acceltotal= []
    disp1 = {}
    disp2 = {}
    # timestart= []
    
    for i in range(len(testnumber)):
        datacm['%d' % (i+1)]= np.array(pd.read_csv('Data_Center_Mass\\accelAndDispl_15Hz\\%d_AccelAndDispl_15Hz.txt' % testnumber[i], header= None, delim_whitespace= True))
        accelcm['%d' % (i+1)]= np.array(datacm['%d' % (i+1)][:,4])[indexesb[i]:indexesf[i]+1]
        timecm['%d' % (i+1)]= np.array(datacm['%d' % (i+1)][:,0])[indexesb[i]:indexesf[i]+1]
        timecmc['%d' % (i+1)]= timecm['%d' % (i+1)]-timecm['%d' % (i+1)][0]
        disp1['%d' % (i+1)] = np.array(datacm['%d' % (i+1)][:,5])[indexesb[i]:indexesf[i]+1]
        disp2['%d' % (i+1)] = np.array(datacm['%d' % (i+1)][:,6])[indexesb[i]:indexesf[i]+1]
        
        
    for i in range(len(testnumber)):
        dm = []
     
        for j in range(len(disp1['%d' % (i+1)])):
            d1 = disp1['%d' % (i+1)][j]
            d2 = disp2['%d' % (i+1)][j]
            dm.append((d1+d2)/2)
            
        dispdacm['%d' % (i+1)] = np.array(dm)*adjfactors[i]
    
    
    data['1']= np.array(pd.read_csv('13_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    data['2']= np.array(pd.read_csv('14_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    data['3']= np.array(pd.read_csv('15_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    data['4']= np.array(pd.read_csv('16_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    data['5']= np.array(pd.read_csv('17_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    data['6']= np.array(pd.read_csv('18_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    data['7']= np.array(pd.read_csv('19_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    data['8']= np.array(pd.read_csv('20_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    data['9']= np.array(pd.read_csv('24_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    data['10']= np.array(pd.read_csv('25_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    data['11']= np.array(pd.read_csv('26_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    data['12']= np.array(pd.read_csv('27_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    data['13']= np.array(pd.read_csv('28_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    data['14']= np.array(pd.read_csv('40_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    data['15']= np.array(pd.read_csv('41_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    
    
    for i in range(15):
        accelbase['%d' % (i+1)]= np.array(data['%d' % (i+1)][:,1])
        timearray['%d' % (i+1)]= data['%d' % (i+1)][:,0]-data['%d' % (i+1)][0,0]
        energyarraycm['%d' % (i+1)]= [np.trapz(accelcm['%d' % (i+1)][:j]*m*(-1),dispdacm['%d' % (i+1)][:j]/1000) for j in range(len(accelcm['%d' % (i+1)]))]
        
    cc= 0
    
    for i in range(len(testnumber)):
        for j in range(len(timearray['%d' % (i+1)])):
            if i==0:
                timetotal.append(timearray['%d' % (i+1)][j])
            else:
                timetotal.append(cc+timearray['%d' % (i+1)][j])
            acceltotal.append(accelbase['%d' % (i+1)][j])
        cc= timetotal[-1]
        if i!=len(testnumber)-1:
            for j in range(int(round(15/timearray['%d' % (i+1)][1]))):
                if i==0:
                    timetotal.append(timearray['%d' % (i+1)][-1]+(j+1)*timearray['%d' % (i+1)][1])
                else:
                    timetotal.append(cc+(j+1)*timearray['%d' % (i+1)][1])
                acceltotal.append(0)
            cc= timetotal[-1]
    
    #### SDOF Analysis
    
    # for i in range(15):
    np.savetxt('TH/Time_Total.txt', timetotal)
    np.savetxt('TH/Acceleration_Total.txt', acceltotal)
    
    recordsa= 'TH/Acceleration_Total.txt'
    recordst= 'TH/Time_Total.txt'
    
    dts= 0.001953
    mass= m/1e6
    
    dur= timetotal[-1]
    
    fh= open('FMdata.tcl', 'w')
    
    fh.write(
        'set FMfilet "%s";\nset FMfilea "%s";\nset dtg %f;\nset durs %f;'
        '\nset M %f;'
        '\nset s1p %f;\nset e1p %f;\nset s2p %f;\nset e2p %f;\nset s3p %f;'
        '\nset e3p %f;\nset pinchX %f;\nset pinchY %f;\nset damage1 %f;'
        '\nset damage2 %f;\nset beta %f;' 
        '\nset s1n %f;\nset e1n %f;\nset s2n %f;\nset e2n %f;\nset s3n %f;\nset e3n %f;'
        '\nset width %d;\nset height %d;'
        '\nset k1 %f;\nset k2 %f;\nset sigAct %f;\nset betaf %f;\nset epsSlip %f;\nset epsBear %f;\nset rBear %f;' % (recordst, recordsa, dts, dur, mass, s1p, e1p, s2p, e2p, s3p, e3p, par[1], par[2], damage1, damage2, beta, s1n, e1n, s2n, e2n, s3n, e3n, width, height, k1, k2, sigAct, par[0], epsSlip, epsBear, rBear))
        # % (recordst, recordsa, dts, dur, mass, s1p, e1p, s2p, e2p, s3p, e3p, pinchX, pinchY, damage1, damage2, beta, s1n, e1n, s2n, e2n, s3n, e3n, width, height))
    fh.close()
    subprocess.call('OpenSees MDOF_CalibrationMH_Grav.tcl', shell=True)
    
    ####
    
    sdofd= {}
    sdofa= {}
    sdoft= {}
    sdoff= {}
    sdofe= {}
    
    dtt = 0.001953
    
    indexes_Beg = [0,int(42/dtt),int(85/dtt),int(127/dtt),int(170/dtt),int(212/dtt),int(255/dtt),int(298/dtt),int(342/dtt),int(387/dtt),int(431/dtt),int(475/dtt),int(518/dtt),int(562/dtt),int(606/dtt)]
    indexes_End = [int(28/dtt),int(71/dtt),int(113/dtt),int(156/dtt),int(198/dtt),int(241/dtt),int(284/dtt),int(328/dtt),int(372/dtt),int(416/dtt),int(460/dtt),int(503/dtt),int(547/dtt),int(591/dtt),324987]
    timei= [0,42,85,127,170,212,255,298,342,387,431,475,518,562,606]
    
    for i in range(len(testnumber)):
        sdofd['%d' % (i+1)]= np.array(pd.read_csv('ResultsAtt/DispT.out', header= None, delim_whitespace= True)[1][indexes_Beg[i]:indexes_End[i]])*0.001
        sdofa['%d' % (i+1)]= np.array(pd.read_csv('ResultsAtt/AccelT.out', header= None, delim_whitespace= True)[1][indexes_Beg[i]:indexes_End[i]])
        sdoft['%d' % (i+1)]= np.array(pd.read_csv('ResultsAtt/AccelT.out', header= None, delim_whitespace= True)[0][indexes_Beg[i]:indexes_End[i]])-timei[i]
        sdoff['%d' % (i+1)]= sdofa['%d' % (i+1)]*0.001*m*-1
        sdofe['%d' % (i+1)]= [np.trapz(sdoff['%d' % (i+1)][:j],sdofd['%d' % (i+1)][:j]) for j in range(len(sdoff['%d' % (i+1)]))]
    
    ErrorEn = 0.0
    ErrorDisp = 0.0
    ErrorFor = 0.0
    
    maxE = np.zeros(15)
    maxD = np.zeros(15)
    maxF = np.zeros(15)
    
    for i in range(1,16):
       maxE[i-1] = max(np.absolute(energyarraycm['%d' % i]))
       maxD[i-1] = max(np.absolute(dispdacm['%d' % i]))
       maxF[i-1] = max(np.absolute(accelcm['%d' % i]))
     
    #Enmax = max(maxE)
    #Dispmax = max(maxD)
    #Formax = max(maxF)
            
            
    for i in range(1,16):
        Net = len(energyarraycm['%d' % i])
        Nen = len(sdofe['%d' % i])
        Ne = min([Net,Nen])
        for j in range(Ne):
            ErrorEn = (energyarraycm['%d' % i][j]/maxE[i-1]-sdofe['%d' % i][j]/maxE[i-1])**2
            ErrorDisp = (dispdacm['%d' % i][j]/maxD[i-1]-sdofd['%d' % i][j]*1000/maxD[i-1])**2
            ErrorFor = (accelcm['%d' % i][j]*m/maxF[i-1]-sdoff['%d' % i][j]/maxF[i-1])**2
            
    wDisp = 0.5
    wFor = 0.1
    wEn = 0.4
    
    TotError = wDisp*ErrorDisp+wFor*ErrorFor+wEn*ErrorEn
    
    return TotError

HystModelFit(params)
    
res = minimize(HystModelFit,params,method = 'Nelder-Mead',options={'disp': True})  

res.x
